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Abstract 

Systematic simulations are carried out based on the model of fluidized beds 

proposed by the present authors [K. Ichiki and H.Hayakawa, Phys. Rev. E 52, 

658 (1995)]. From our simulation, we confirm that fluidization is a continuous 

transition. We also confirm the existence of two types of fluidized phases, the 

channeling phase and the bubbling phase. We find the close relations between 

the averaged behaviors in fluidized beds and quasi equilibrium states in dense 

liquids. In fluidized beds, (i) the flow rate plays the role of the effective 

temperature, and (ii) the existence of a kind of the fluctuation-dissipation 

relation is suggested. 
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I. INTRODUCTION 



Recently granular materials have been studied extensively from both experimental and 
theoretical point of views in the context of the nonequilibrium statistical physics Since 
the granular materials are dissipative, energy injections are necessary to preserve steady 
states. Many of recent studies for granular materials are focused on the behavior of systems 
excited by mechanical activations such as vibration or rotation of vessels. On the other 
hand, the researches on fluidized beds, where systems are excited by the fluid flow, are not 
relatively advanced in spite of their variety of dynamical behaviors H||. 

Fluidized beds are widely used in chemical industries since early 19 's, and they have 
been studied from technological point of views. Fluidized beds consist of granular particles 
confined in a tall chamber with distributor for the fluid flow at the bottom. In experiments, 
energy injection to the system is controlled by the flow rate of fluid. At low flow rate, 
the system is in the fixed phase where particles rest on the bottom. When the flow rate 
exceeds the critical value, particles start moving. This state is called the fluidized phase, 
which contains sub phases, for instance, the homogeneous phase, the bubbling phase, the 
channeling phase, etc. 

There are many models to describe fluidized beds, which can be classified in two cat- 
egories, two-fluid models and particle-dynamics models. In the two-fluid models, particles 
are treated as a fluid flBHlOll. These models have benefit on analytical treatments and to 



generalize their discussion to other systems pT| . However, their bases, such as constitution 
equations for the particle-phase pressure and the stress tensor, have not been established. On 
the other hand, the particle-dynamics models describe the direct motion of particles. There 
are various models, which are kinetic theories ||, the discrete element method |T^-[14|, etc. 
However, these models cannot be the basis for the two-fluid models. The main problem of 
them is that hydrodynamic interactions among particles are over simplified. For instance, 
the boundary condition between particles and fluid is not satisfied in the scale of particles, 
and the fluid equation is calculated under inviscid limit. 
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Recently the present authors have proposed a numerical model of fluidized beds, where 



hydrodynamic interaction among particles is calculated with reliable accuracy |I3 , P!6| . In this 
paper, we will present the results of our systematic simulations and behaviors on statistical 
quantities obtained from the simulations. 

The contents of this paper are as follows. In Sec. II, we review the method of our 
simulation. We show the results in Sec. Ill, where we observe the transition of fluidization 
and the existence of two fluidized phases. We also discuss statistical quantities, which are 
analogous to equilibrium correspondences. In Sec. IV we give an interpretation of the 
averaged quantities by the hole theory for simple liquids [|17]>II1] • In Sec. V, we conclude our 
results. In Appendices, we summarize the method of our simulation and discuss theoretical 
difficulties in the modeling of fluidized beds. 



II. SIMULATION METHOD 

In this section, we briefly explain our model and how to simulate the dynamics of granular 
particles in fluid flows. The detail explanation of our model can be seen in Ref. [[H|[nj and 
Appendices. For simplicity, we only consider the cases of monodispersed spherical particles. 
We assume the following equation of motion for the particles 

St^ = -U + V + F c , (2.1) 
dt 

where St is the effective Stokes number, F c represents hard-core collision among particles, 
U is the velocity of particles. V is the terminal velocity determined by 

V - u°° = -FT 1 ■ E z , (2.2) 

where u°° is the flow rate of induced fluid which is equal to the superficial velocity conven- 
tionally used for the fluidized beds. R is the resistance matrix representing the hydrody- 
namic interaction among particles calculated by the method of the Stokesian dynamics ||19[ . 
where periodic boundary condition is adopted as the effect of chamber (see Appendix |A2j ). 
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Hard-core collisions are assumed to be elastic and calculated by the momentum exchange 
for contacting particles in simulations (see Appendix [A 1| ). The bold-face letters without 
superscripts represents vectors in 3iV-dimension, where N is the number of particles in the 
unit cell of periodic boundary condition. For example, the velocity U has the following 
components 



U 



U(2) 
IJ(N) 



(2.3) 



where the bold-face letters with superscripts represent vectors in 3-dimension. In this paper, 
we use dimensionless quantities with the aid of the particle radius a and the sedimentation 
velocity of a single particle in a viscous fluid Uq = mg/Qnpa, where /i is the viscosity 
of the fluid, m is the mass of the particles and g = g(p p — Pf)/p P with the gravitational 
acceleration g, and the densities of the particle p p and the fluid pf. Equation ( |2.1|) represents 
the relaxation process of U to V with the time-scale St. 

In Eq. Q2.1|) , there are two control parameters, the effective Stokes number St and the 
flow rate u°°. For the parameters related to the system size, we choose the number of mobile 
particles A^m = 256, that of fixed particles N-p = 10 and the size of the unit cell in periodic 
boundary condition (L x , L y , L z ) = (34,2,100). In this situation, particles are confined in 
the vertical plane, while hydrodynamic interactions are considered in 3-dimensional space. 
We adopt the fixed phase as initial conditions of our simulations, which is constructed from 
simulations with m°° = 0. The choice of the system size and these artifact situations come 
from the limitation of computer resources. We have checked that statistical quantities seems 
to be insensitive to the choice of L z within the range of 50 < L z < 100 and the choice of the 
initial conditions is not relevant from the comparison of results with other initial conditions. 
We have also confirmed that qualitatively similar behaviors to this situation are observed in 
3-dimensional simulations and in the case of A^m = 133. 
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III. SIMULATION RESULTS 



A. General behavior 

In this section, we present the results of our simulations in details. We perform simu- 
lations at the points in the parameter space (Fig. |l|) within the range of 0.05 < u°° < 0.8 
and 0.1 < St < 100, where we observe fixed, bubbling and channeling phases. Transitions 
among these phases will be discussed below. 

In the fixed phase at low flow rate, particles are rest at the bottom. At the critical flow 
rate u°° = u c , the particles begin to be fluidized. It seems that the transition between the 
fixed phase and the fluidized phase is independent of St. We observe two fluidized phases. 
One is the channeling phase observed for small St where we can see a channel or a path of 
fluid flow. Another is the bubbling phase observed for large St where bubbles raise through 
the particle beds. Typical snapshots of them are shown in Figs. and [3]. We show the area 
of channel-bubble transition observed in our simulations as the transition area in Fig. [l], 
where we will see that statistical quantities qualitatively change their behaviors. 

To characterize the transition of fluidization quantitatively, we calculate the kinetic en- 
ergy per particle E(t) defined by 

I N M 

E(t) = ^Y.\V a (t)\ 2 - (3.1) 

A typical behavior of E(t) in the bubbling phase is shown in Fig. f|. We observe regular 
behavior after the minimum point following the first peak of E(t). In the channeling phase, 
we also observe qualitatively similar behavior of E(t) to Fig. |], though the period of peaks is 
smaller and each peak is not distinguishable. We now introduce the average of E(t) defined 
by 

where AT is the period of regular behavior in E(t). From Fig. |5|, we observe a continuous 
transition at u°° = u c and the linear behavior of E(u°°) in the fluidized phase (u°° > u c ). 
This behavior can be well fitted by 
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(u°° < u c ) 

E(u°°) = { , (3.3) 

A E (u°°-u c ) (u°°>u c ) 



where A E and u c are the fitting parameters which depend on St. Equation ( |3.3j ) defines u c 
shown in Fig. [l]. 

Now we show that the transition of fluidization can be understood as the process gen- 
erating the free volume around the particles in the fixed phase. It is useful to remember 
that our model is Galilei invariant, that is, the system with the fixed particles of Uf = 
under the flow rate u°° is equivalent to that with the fixed particles of Uf = — m°°E 2 under 
the flow rate 0. Let us consider the process under the latter situation. First we define 
which is the falling velocity of the mobile particles in the fixed phase without the support 
of fixed particles. If the flow rate u°° is smaller than {/fail, the mobile particles cannot pass 
over the fixed particles moving downward with — u°°. Therefore the mobile particles hold 
on the fixed particles and the gap between the mobile particles and the fixed particles is not 
generated. While the flow rate u°° is larger than {/fail, the mobile particles apart from the 
fixed particles and then the gap between them is generated. The gap causes the Rayleigh- 
Taylor instability observed when a heavy fluid exists above a light fluid [[HJ • Then the gap 
may grow into a bubble and propagate upward through the particles, or may construct a 
channel. From this discussion, the critical flow rate u c is determined by the falling velocity 
{7f a ii. This suggests that u c is independent of St because the falling velocity {/f a n can be 
evaluated as the sedimentation rate of suspensions [21 . 

Next we discuss the channel-bubble transition. In view of Figs. and |3|, it is hard to 
distinguish the channeling phase from the bubbling phase. At first, we show the variance 
Vh defined by 

Vh = ^T J AT dt (H(t) ~ H)\ (3.4) 
where H and H are the height of the center of mass and its average defined by 

i N M 
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and 

B = L* "<■<>■ (3 - 6) 

We expect that Vh is small in the channeling phase and large in the bubbling phase. Figure 
|6] shows the corresponding behavior for u°° = 0.3 and the transition is observed around 
St — 5. For other cases, the channel-bubble transitions are observed in the area shown in 

Fig. a 

We also observe some qualitative changes in the physical quantities corresponding to the 
channel-bubble transitions. Here we discuss the behavior of E(St). From Fig. [5], we see that 
E increases with St in the channeling phase, and E decreases with St in the bubbling phase. 
The behavior in the channeling phase can be understood by the following reason. At St = 0, 
we observe a steady channel, no relative motions among particles. When St increases, the 
particles on the channel can be fluidized. This is because St is the relaxation time of the 
particle velocity U to their terminal velocity V, which is determined under the case of 
St = 0, and the non-relaxed particles cause the collapse of the channel. Therefore, E(St) 
increases with St in the channeling phase. On the other hand, we observe that E decreases 
with St in the bubbling phase. This can be understood as follows. From the previous paper 
15fl or in Fig. |], we have seen that the relative motion of particles or the kinetic energy is 



generated by bubbles which have the difference of the local volume fraction in the system. 
In fact, bubbles are always accompanied by the convective motion of particles. Because St 
is the response time of the velocity for the change of configurations, bubbles become obscure 
as St increases. In fact, we see a sharp bubble and the definite convective motion of particles 
for St = 10 (Fig. H), while we see a relatively obscure bubble and the weaker convection for 
St = 100 (Fig. |9|). Therefore E(St) in the bubbling phase decreases with St. As a result, 
E(St) has a peak around the transition point. This reflects on the St dependence of various 
properties such as fj, e (St) shown in Sec. [Ill B] and H(St) discussed in Sec. [TV[ 
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B. The analogy to equilibrium systems 



In this section, we demonstrate the existence of surprising correspondences in statistical 
quantities between fluidized beds which is in highly nonequilibrium states and quasi equi- 
librium systems. The result of E(u°°) in Fig. |5| suggests that the flow rate u°° behaves as 
the effective temperature of the environment such as that of the heat bath for equilibrium 
systems. Therefore, the critical flow rate u c may correspond to the critical temperature and 
E may be the order parameter of the ffuidization. In the following, we will interpret the 
results of our simulations using this effective temperature u°°. 

In nearly equilibrium systems at temperature T, the resistance £ of a tracer particle is 
given by 

C = -g-, (3.7) 



where D is the diffusion constant and k is the Boltzmann constant p2] . This is the Einstein 
relation and is the simplest form of the fluctuation-dissipation theorem which relates the 
correlation functions in the equilibrium state and the transport coefficients. 
Similarly, we introduce for our simulations the effective viscosity \x e as 

= 7p (3.8) 

Up 

where D p is the diffusion constant of our simulations defined by 

1 1 Nyi 

A> = ^ElU> = 0)| 2 , (3.9) 

1 7V M Q=l 

where XJ a (uj) is the Fourier transform of XJ a (t) calculated by the standard FFT. Equation 
( f378|) can be regarded as an extension of Eq. (|3.7|) with the aid of the effective temperature 
u°°, because the viscosity is usually proportional to the resistance. The observed viscosity 
\i e defined by ( ft.7| ) is shown in Figs. |H] and [11]. 

If the definition of the viscosity ( |3.8| ) is self-consistent, the Einstein relation or the 
fluctuation-dissipation relation in fluidized beds may be valid with the replacement of kT 
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by u°°. This statement is interesting because the system is in a highly nonequilibrium state 
and there is no reason of the existence of the fluctuation-dissipation relation in the sense of 
linear nonequilibrium statistical mechanics. 

We find that the flow-rate dependence of the viscosity fj, e (u°°) obtained by Eq. fl3.8p 
obeys the Arrhenius function 

He(u°°) oc e £/M °°, (3.10) 



where e is a fitting parameter. The fitting by Eq. ( p. 10 ) with e = 0.113 ±0.017 is also shown 



in Fig. 0. We compare our result of /i e with experimental result in fluidized beds [23]. In 
the experiment of fluidized beds, the shear viscosity measured by the modified Stormer 
viscometer also obeys the Arrhenius function of Eq. fl3.10 ) [23]. Therefore, our result from 



Eq. ( j3.8|) is consistent with the experiment. This behavior which can be understood by a 
dense liquid theory in part will be explained in Sec. 

The connection between the non-Gaussian property and the dissipation in the system 
have been discussed for granular materials [f2^,^,|l^]. Here we check the non- Gaussian 
property in the velocity distribution functions P(U X ) obtained from our simulations. To 
characterize the non-Gaussian property of the velocity distribution, we calculate the 4th 
cumulant C 4 defined by 

Ci(U s ) = (U x 4 ) - 3(U X 2 ) 2 - 4(U X )(U X 3 ) 

+ 12(U X ) 2 (U X 2 }-Q(U X } 4 . (3.11) 



The resultant behaviors of C4 are shown in Figs. [12] and [13], where they are scaled by the 
square of variance (or 2nd cumulant) defined by 

C 2 (U X ) = (U x 2 ) - (U x ) 2 . (3.12) 

This non-Gaussian parameter C^jiC-i) 2 is zero for the Gaussian distribution and 3 for the 
exponential distribution. These behaviors of C^j{C%f in Figs. O and O are similar to 



those of the effective viscosity fi e in Figs. |TU] and [II]. In fact, we can also fit C^jiC?) 2 by 
the Arrhenius function 



C 4 /(C 2 ) 2 oc exp^'/w 00 ), (3.13) 

as shown in Fig. O with e' = 0.175 ± 0.045. It will be an interesting subject that we will 
study a quantitative relation between (13.10Q and ( |3.13|) . 



IV. DISCUSSIONS 

In this section, we demonstrate that the qualitative understanding of the results of our 
simulation is possible with the aid of the hole theory applied to simple liquids. 

First we discuss the flow rate dependence of the viscosity /j, e (u°°). The hole theory, which 
is used for the behavior of simple liquids, is based on the following picture. A molecule in a 
liquid can move when a free volume or a hole is generated around it. From this picture, the 
empirical relation of the viscosity /!/ |Oj is derived as 



H(T) oc exp , (4.1) 

where E\ is the activation energy to create the free volume. If we regard u°° as the temper- 
ature, Eq. ( |3.10| ) is identical to Eq. ( |4.1| ). This suggests that the averaged behavior in the 



fluidized beds may be understood by the hole theory in simple liquids p6| - 

Next we discuss St dependence of the height of center of mass H(St) defined by ( [3.61) . 



In Fig. |lj we show a typical behavior of H(St) at u°° = 0.3. We see the qualitative change 
of behavior in the transition area in Fig. p], where H is almost constant in the channeling 
phase, and H increases logarithmically in the bubbling phase. Since the change of behavior 
with St in the channeling phase is only how the channel collapses, H(St) is expected to be 
independent of St. While the behavior in the bubbling phase is interesting. We can fit the 
data in Fig. |TJ as 

H(St) = C H \og(St) + D H , (4.2) 



where the fitting parameters Ch and Du depend on u°° . Equation ( f4.2|) can be rewritten as 



St = exp [SjlE^j . (4.3) 
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Here H — Dh can be understood as the volume expansion AV, since Dh is the height at 



St = 1. We show u°° dependence of Ch in Fig. |TJ|. From this figure, Cjj is approximately 
proportional to the effective temperature u°° . Since St is the characteristic time r of the 
system, we rewrite Eq. ( |4.3| ) as 

T K ex P > ( 4 - 4 ) 



where F is a constant. Equation (4.4) is consistent with Eq.( fOj) under the reasonable 



assumption where the viscosity is characterized by the time r to generate the free volume, 
and the activation energy to generate the free volume or expansion is proportional to AV. 
Before closing this section, we give some remarks. Although the transition of fluidization 



in experiments seems to be the discontinuous phase transition , our simulations suggests a 
continuous phase transition. Also it is an open problem that at present we cannot reproduce 
homogeneous phase in our simulation. For these problems, we need to examine carefully 
the difference between the experiments and the simulations, and we must investigate the 
behavior near the critical flow rate u c in detail, because the discontinuous phase transition 
and the homogeneous phase are observed there in experiments. 



V. CONCLUSIONS 

In this paper, we have carried out systematic simulations with the change of two control 
parameters, the flow rate of the fluid u°° and the effective Stokes number St. When the flow 
rate u°° is small, particles rest in the fixed phase. Above the critical flow rate u c , particles 
are fluidized. The critical value u c is independent of St. We have found two fluidized phases, 
the channeling phase and the bubbling phase, where the former changes to the latter as St 
increases. 

We have found that the flow rate u°° plays the role of the effective temperature. In 
terms of the effective temperature u°°, we have defined the effective viscosity fi e with the 
aid of the Einstein relation. The flow-rate dependence of the viscosity fi e is similar to that 
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in the experiments in real fTuidized beds. We also find that the viscosity fj, e (u°°, St) can be 
an index of the non-Gaussian property in velocity distribution of particles. This property 
is consistent with the behavior on granular materials or the system of inelastic particles. 
Qualitative behavior of fluidized beds such as fj, e (u°°) and H(St) can be understood by 
means of the hole theory which has been used for simple liquids. 
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APPENDIX A: THE IDEAL MODEL OF FLUIDIZED BEDS 



Here we review our model of fluidized beds, presented in the previous paper | 15| , where 
only essences are extracted from real systems and all other irrelevant mechanisms are ne- 
glected. 



1. Equation of motion 

We construct the model by only four mechanisms, which are the inertial effect of the 
particles, hydrodynamic interaction through the fluid, the gravitational force and the contact 
force. Therefore the equation of motion can be written as 

5t ^ = F / + F g + F c , (Al) 

where Fj,F 9 ,F c are the force from the fluid, the gravitational force and the contact force 
respectively. St is the bare Stokes number defined by 
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St a = (A2) 

Particles are assumed to be monodisperse and hard-core spheres and rotational motions and 
torques acting on particles are neglected. 

The gravitational force F g can be written as 

F 9 = -E z , (A3) 

where E z is the unit vector directed to the z axis by the notation of Eq. (|2.3|) . It is assumed 
that the direction of the gravity is —z. 

It is assumed that F c is the impulse by prefect elastic collisions. Because collisions are 
inelastic in real systems, this assumption of elastic collision means our standpoint of model- 
ing that the essential mechanism of fluidized beds is not the inelasticity in collisions but the 
hydrodynamic interaction. Even in this case, the model is dissipative, because the hydro- 
dynamics interaction is nothing less than the friction. In our simulation, we represent the 
collision by the momentum exchange at contact in stead of the contact force F c . Therefore 
we do not write F c explicitly in the equations in the following discussion. 



2. Hydrodynamic interaction 

In our model, hydrodynamic interaction among particles through the fluid is considered 
under the Stokes approximation where the viscous effect of the fluid dominates the inertia 
of the fluid. The reason to adopt the Stokes approximation is as follows: The hydrodynamic 
interaction is the friction between the particles and the fluid and the friction is originated 
by the viscosity of the fluid. 

In the Stokes approximation, the force acting on the particles from the fluid F/ and the 
velocity of the particles U are related by the resistance matrix R as 

-F f = R- (U -u°°), (A4) 

where u°° is the velocity of the fluid without particles. R contains all information about the 
interaction and depend only on the configuration of particles. 
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To calculate R, we adopt the method in the Stokesian dynamics, which is developed by 



J.F.Brady and his collaborators for dense colloidal particles ||30|j31| , p!9|| . In the Stokesian dy- 
namics, R is constructed from the two contributions in limited cases, which are the mobility 
matrix in dilute limit M and the exact resistance matrix in two-body problem R2B [|32 
as follows, 



R=(M )- x + R 



lub 



(A5) 



♦lub 



where i? is constructed by the pairwise-additive manner from the two-body lubrication 



+lub 



matrix i? 2B which is defined by 



♦lub 



R2B - R2B - (M 2B ) . 



(A6) 



M 2B is the two-body mobility matrix in dilute limit. In general M is formulated by the 
multipole expansion P0f . In the model of fluidized bed, however, we need to introduce the 
effect of the chamber. The chamber in the real fluidized beds has two contributions which 
are to bound the fluid by the vertical wall and to support the particles by the bottom. 

We introduce the contribution of bounding the fluid in terms of the periodic boundary 
condition. Because M has the long-range interaction, we use the Ewald summation tech- 
nique f33|. We can construct the resistance matrix under the periodic boundary condition 

«-» 00 

also by Eq. ( |A5| ) only replacing M to the Ewald summed tensor for N particles in the 

unit cell |H] . 



On the other hand, the contribution of the bottom supporting the particles is introduced 
by the particles fixed in space. In this case, we get the force from the fluid to the mobile 
particles Ff in Eq. ( |A1| ) as follows, 

— Ff — Rmm • (Um — u°°) — R M f ■ u°°, (A7) 

where the subscripts "M" and "F" represent mobile and fixed particles respectively. The 
complete form of the resistance relation is 

TT,„-ii°° 

(A8) 



F/ 




Rmm 


<— > 

Rmf 






F F 




Rfm 


Rff 




- u°° 
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where Ff is the force acting on the fixed particles. 

For simplicity we only discuss the case without fixed particles in the following. However 
we can get the correct forms only the replacement of Ff of (|A4| ) by ( |A7| ) . 



3. Effective inertia 

The inertia of particles causes the relaxation process on the velocity from the initial value 
to the optimal value where the inertia corresponds to the relaxation time. The optimal value 
is usually called the terminal velocity determined on the steady state. For the fluidized beds, 
the terminal velocity V is determined by Eq. (|A1|) with dXJ/dt = and get Eq. (|2.2|) . This 
is the case that all forces acting on particles balance. With this terminal velocity we can 
write the relaxation process as Eq. (|2.1|) . 



Equation ( |2.1| ) is the same as used in the previous paper JTSJ. Although we had argued 
that this equation might be justified in some approximation, we would state here that fl2.ip 
contains all of essential processes in fluidized beds. On this point, we will discuss in Appendix 



APPENDIX B: THE INERTIAL EFFECT OF PARTICLES 

Here we discuss the difficulties arising when we use Eq. ( |A1|) with Ff in the Stokes 
approximation (|A|). 

If we write Eq. ( |Al|) with Ff, we get 

StoR' 1 ■ j\J = -U + V. (Bl) 



From the simulation of Eq. (|B1|) , we observe no collision between particles even in the case 
with large St . Particles form a cluster and relative motions among them almost disappear. 
This situation may be understood by the following model, 

a »f = -^ K (B2) 
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Here we extract the radial component of the motion between two particles, r denotes the 
separation between the centers of the pair and U denotes the relative velocity. The resistance 
l/(r — 2a) reflects the lubrication effect which diverge at contact (r — > 2a). In this model, 
we need initially the infinite energy to approach the contact point (r = 2a) even in the case 
of large St^. Thus the model (|B2[) cannot contains any collision. 



Our model fl2.1|) can be understood as the renormalization of the singularity in the 
lubrication, because we get Eq. ( |2.1|) by multiplying the singularity R on the inertial term of 
( pip . (Here we note that R is dimensionless because it scaled by 6717/a.) Our model behave 
reasonably like the real fluidized beds where collisions are occurred so frequently. This 
suggests that the singularity of the lubrication must be prevented from some mechanisms in 
the real systems. 

We can imagine several possible mechanisms preventing the singularity. For example, if 
there are some dimples on the surface of the particles, they collide before the mean surfaces 
contact. From another point of view we can also say that the continuous description of the 
fluid in the gap between the particles breakdowns when particles approach closely and the 
gap becomes comparable to the mean-free path of the molecules of the fluid 



Recently a model in this context has been presented |35|| . They introduce a cut-off length, 
which may correspond to the height of the dimples on the surface or the mean-free path 
of the fluid molecules. If particles approach with each other within the cut-off length, the 
gap between them is assumed to be the cut-off in calculation of hydrodynamic interaction. 
The result of their simulation is suggestive even though the situation, which is the behavior 
in the shearing flow without the gravity, is different. Their results are characterized by the 
parameter St s /(R C ), where St s = rwy/Gn^a is the Stokes number in shear flow with the 
shear rate 7 and (R c ) is the averaged resistance of a tracer with the same volume fraction 
and the cut-off length. 

From the above discussion, we get the meaning of the effective Stokes number St as 

St = St /(R c ), (B3) 
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which depends on the cut-off length of the real systems. We need more delicate investigations 
for the dependence of the cut-off length or the mechanism preventing the singularity in the 
lubrication. 
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FIGURES 

FIG. 1. Simulations performed in the parameter space (u°°,St) with Nm = 256, Np = 10 and 
(L x , L y , L z ) = (34,2, 100). In this figure we show the observed phases, fixed phase (•), channeling 
phase (□), bubbling phase (o) and transition phase (A). We also show the transition line of 
fluidization u c with solid line and the area of the channel-bubble transition with dotted line, which 
are discussed in the text. 

FIG. 2. Typical snapshots of channeling phase with St = 0.5, u°° = 0.15. The time proceeds 
from left to right with the interval 20 dimensionless time. 

FIG. 3. Typical snapshots of bubbling phase with St = 10.0, u°° = 0.15. The time proceeds 
from left to right with the interval 20 dimensionless time. 

FIG. 4. A typical behavior of E(t) in the fluidized phase. Parameters used are 
St = 20, u°° = 0.3. The oscillation corresponds to the generation of bubbles. 1 step means 1 
dimensionless time. 

FIG. 5. The flow rate dependence of the averaged energy E(u°°) with St = 10. Error bars are 
their standard deviation. We can see that the transition of fluidization at u°° = u c and the linear 
behavior of E for the flow rate u°° in the fluidized phase (u°° > u c ). 

FIG. 6. St dependence of the variance Vh with u°° = 0.3. We can see the transition around 
St = 5. At the channeling phase the variance is small, while at the bubbling phase the variance 
become large. 

FIG. 7. St dependence of the averaged energy E(St) with u°° = 0.3. We can observe a peak 
around the channel-bubble transition shown in Fig. [l]. 

FIG. 8. The convective motion of particles in the bubble with u°° = 0.3, St = 10. Here two 
periodic images are shown. We can observe the sharp edge of the bubble and the definite convection. 
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FIG. 9. The convective motion of particles in the bubble with u°° = 0.3, St = 100. Here two 
periodic images are shown. Comparing to the Fig. ||, we can observe broader edge of the bubble 
and weaker convection, where the scale of the velocities are the same. 

FIG. 10. Effective viscosity fj, e (u°°). This result is calculated on the simulation with St = 10.0. 



The fitting by the Arrhenius type function ( 3.10 ) with e = 0.113 ± 0.017 is also shown. 



FIG. 11. Effective viscosity fj, e (St). This result is calculated on the simulation with u°° = 0.3. 

FIG. 12. The non-Gaussian parameter (^/(C^) 2 with St = 10 for the change of u°° . We also 
show the fitting by exp(e'/u°°) with e' = 0.175 ± 0.045. 

FIG. 13. The non-Gaussian parameter C^jiC-i) 2, with u°° = 0.3 for the change of St. 



FIG. 14. St dependence of H with u°° = 0.3. For the bubbling phase, we fit by Eq.(4.2) with 
C H = 3.087 ± 0.055, D H = 16.51 ± 0.14. 

FIG. 15. u°° dependence of the fitting parameter Ch- We can see Ch is linear for u°°. 
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